function [LE1, LE2, xStore] = LE(f, d, x0, phi0, N_iter, Transient, KeepNum)
% LE_improved_COGM 计算 CO-GM 映射的两个李雅普诺夫指数
% 输入:
%   f         - 函数句柄, f(d, x, phi, i) 返回 [x_next, phi_next]
%               其中 i 为全局迭代计数
%   d         - 映射参数（耦合参数，CO-GM 中 d 对应 f 的第一个参数）
%   x0,phi0   - 初始条件
%   N_iter    - 总迭代步数
%   Transient - 舞台瞬态步数（不计入 LE 累计）
%   KeepNum   - 用于记录分岔图数据的最后步数
% 输出:
%   LE1, LE2  - 两个李雅普诺夫指数（保证 LE1>=LE2）
%   xStore    - 记录最后 KeepNum 步的 x 值

    % 初始化状态与扰动向量
    x   = x0;
    phi = phi0;
    v1 = [1; 0];
    v2 = [0; 1];
    sumLog1 = 0;
    sumLog2 = 0;
    
    % transient 阶段：更新状态，不累计伸长因子
    for i = 1:Transient
        [x, phi] = f(d, x, phi, i);
    end
    
    % 用于保存分岔图数据
    xStore = zeros(1, KeepNum);
    idxStore = 1;
    
    % 主循环：全局迭代计数从 Transient+1 到 N_iter
    for i = Transient+1 : N_iter
        % 1. 更新扰动向量
        [v1, v2] = updateVectors(f, d, x, phi, i, v1, v2);
        
        % 2. 计算并累计伸长因子（未归一化前的模长）
        normv1 = norm(v1);
        normv2 = norm(v2);
        sumLog1 = sumLog1 + log(normv1);
        sumLog2 = sumLog2 + log(normv2);
        
        % 3. 对扰动向量归一化，并做 Gram-Schmidt 正交化
        v1 = v1 / (normv1 + 1e-30);
        proj = dot(v2, v1) * v1;
        v2 = v2 - proj;
        normv2 = norm(v2);
        if normv2 ~= 0
            v2 = v2 / (normv2 + 1e-30);
        end
        
        % 4. 更新系统状态
        [x, phi] = f(d, x, phi, i);
        
        % 5. 保存最后 KeepNum 步的 x 值用于分岔图
        if i > (N_iter - KeepNum)
            xStore(idxStore) = x;
            idxStore = idxStore + 1;
        end
    end
    
    N_main = N_iter - Transient;
    LE1 = sumLog1 / N_main;
    LE2 = sumLog2 / N_main;
    
    if LE1 < LE2
        temp = LE1;
        LE1 = LE2;
        LE2 = temp;
    end
end